{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "baba5557",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.stats import norm\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e55631a",
   "metadata": {},
   "source": [
    "# Identification of the parameters $\\{p, \\phi, \\alpha_0, \\alpha_1\\}$\n",
    "\n",
    "The model is\n",
    "$$\n",
    "\\max_{m_i} \\frac{(y_i-pm_i)^{1-\\sigma}}{1-\\sigma} + \\phi (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "The FOC is \n",
    "$$\n",
    "p (y_i-pm_i)^{-\\sigma} = \\phi\\alpha_0 \\exp(-\\alpha_1 - \\alpha_0 m_i)\n",
    "$$\n",
    "This gives a solution $m_i = {\\cal M}(y_i,p|\\sigma, \\phi, \\alpha_0, \\alpha1)$\n",
    "\n",
    "\n",
    "### US\n",
    "For the US, we have $p=1$ and $\\sum_i \\omega_i y_i \\equiv y = 1$. Therefore, the moments restrictions  \n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i}\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "$$\n",
    "\\Psi_3 = \\frac{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\max}))}{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\min}))}\n",
    "$$\n",
    "allows to identify $\\{\\phi,\\alpha_0,\\alpha_1\\}$\n",
    "\n",
    "\n",
    "### Europe\n",
    "\n",
    "For Europe, we keep $\\{\\phi,\\alpha_0,\\alpha_1\\}$ found in the previous step, and the moments restrictions\n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i}\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "allows to identify $\\{p,\\alpha_1\\}$, given that a calibration of $\\sum_i \\omega_i y_i \\equiv y$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e3d32d8",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5f76a0af",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function for finding the decision rule for m\n",
    "\n",
    "def medexp(m, *param_m):\n",
    "    p, y, sigma, alpha0, phi, alpha1 = param_m \n",
    "    return p*( (y - p*max(0,min(m,y/p)))**(-sigma) ) - alpha0*phi*np.exp(-( alpha1 + alpha0*max(0,min(m,y/p)) )) \n",
    "\n",
    "# Function for finding parameters (alpha0, phi, alpha1) using US moments' restrictions\n",
    "def hetero_calib(pd_c, *param_c):\n",
    "    \n",
    "    Psi1, Psi2, Psi3, p, yy, sigma = param_c \n",
    "    \n",
    "    msolv = np.zeros(yy.size)\n",
    "    \n",
    "    for ii in range(yy.size):\n",
    "        param_cm  = (p, yy[ii], sigma, pd_c[0], pd_c[1], pd_c[2])\n",
    "        m0        = .1\n",
    "        msolv[ii] = max(0,fsolve(medexp, m0, args=param_cm))\n",
    "    \n",
    "    return [Psi1 - np.mean(msolv), \n",
    "            Psi2 - np.mean(1-np.exp(-pd_c[2]-pd_c[0]*msolv)),\n",
    "            Psi3 - (1-np.exp(-pd_c[2]-pd_c[0]*msolv[-1]))/(1-np.exp(-pd_c[2]-pd_c[0]*msolv[0])) ] \n",
    "\n",
    "# Function for finding (p, alpha1) using European moments' restrictions\n",
    "def hetero_calib_eu(pd_s, *param_s):\n",
    "    \n",
    "    Psi1, Psi2, yy, sigma, phi, alpha0 = param_s \n",
    "    \n",
    "    msolv = np.zeros(yy.size)\n",
    "    \n",
    "    for ii in range(yy.size):\n",
    "        param_sm  = (pd_s[0], yy[ii], sigma, alpha0, phi, pd_s[1])\n",
    "        m0        = .1\n",
    "        msolv[ii] = max(0,fsolve(medexp, m0, args=param_sm))\n",
    "    \n",
    "    return [Psi1 - (pd_s[0]*np.mean(msolv))/(np.mean(yy)), \n",
    "            Psi2 - np.mean(1-np.exp(-pd_s[1]-alpha0*msolv))] \n",
    "\n",
    "# Function for finding the price, given all other parameters (alpha0, phi, alpha1)\n",
    "def hetero_new_price(p_new, *param_s):\n",
    "    \n",
    "    Psi1, yy, sigma, phi, alpha1, alpha0 = param_s \n",
    "    \n",
    "    msolv = np.zeros(yy.size)\n",
    "    \n",
    "    for ii in range(yy.size):\n",
    "        param_sn  = (p_new, yy[ii], sigma, alpha0, phi, alpha1)\n",
    "        m0        = .05\n",
    "        msolv[ii] = max(0,fsolve(medexp, m0, args=param_sn))\n",
    "    \n",
    "    return Psi1 - (p_new*np.mean(msolv))/(np.mean(yy)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8bf25d4",
   "metadata": {},
   "source": [
    "## Calibrated parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "49f43316",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma  = 2.0\n",
    "Ncount = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf29370a",
   "metadata": {},
   "source": [
    "## Aggregate moments\n",
    "\n",
    "Health inequalities: we use the gradient for the 4th quartile relative to the first. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4acc8877",
   "metadata": {},
   "outputs": [],
   "source": [
    "grad_c  = [1.061, 1.276]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0365aed2",
   "metadata": {},
   "source": [
    "GDP share of health expenditures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ac3df66c",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_c  = [0.0955, 0.1504]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b855728",
   "metadata": {},
   "source": [
    "Income"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "20c3dd9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_c  = [0.78, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa65a7cb",
   "metadata": {},
   "source": [
    "For health status, we use the average health from the ADL definition in Table 4. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "be08d6d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "h_c  = [0.931, 0.892]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a1eab219",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.0955, 0.931 , 1.061 ],\n",
       "       [0.1504, 0.892 , 1.276 ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi = np.zeros((Ncount,3))\n",
    "for ii in range(Ncount):\n",
    "    psi[ii,:] = [s_c[ii],h_c[ii],grad_c[ii]]\n",
    "\n",
    "psi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d01473d0",
   "metadata": {},
   "source": [
    "## Income inequalities: the US\n",
    "\n",
    "Table 5 gives the variance of the income process. If log income has mean 0 and variance equal to $\\sigma^2_{\\epsilon}$, we can solve for the quartiles (we use mid-points of quartiles, 12.5, 37.5, 62.5 and 87.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48a95d99",
   "metadata": {},
   "source": [
    "We take the mean for EU in Table 5. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e4046ef5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.19111728571428574, 0.486429]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sig_eu  = np.mean([0.256367, 0.094643, 0.214538, 0.237439, 0.169834, 0.094643, 0.270357])\n",
    "sig_c   = [sig_eu, 0.486429]\n",
    "sig_c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e5a0143",
   "metadata": {},
   "source": [
    "We find the percentiles, transform in levels and normalize to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a9f5396e",
   "metadata": {},
   "outputs": [],
   "source": [
    "qs   = [0.125,0.375,0.625,0.875]\n",
    "yy_c = np.zeros((Ncount,4))\n",
    "\n",
    "for ii in range(Ncount):\n",
    "    yy    = [norm(0,np.sqrt(sig_c[ii])).ppf(q) for q in qs]\n",
    "    yy    = np.exp(yy)\n",
    "    yy_c[ii,:] = y_c[ii]*yy/np.mean(yy)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9017749",
   "metadata": {},
   "source": [
    "Here are relative incomes. Much more inequality in US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "806c601a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4.975918497985015, 2.7340817164336597)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yy_c[-1,-1]/yy_c[-1,0], yy_c[0,-1]/yy_c[0,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f2a9f4f",
   "metadata": {},
   "source": [
    "## Solution for the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6282c357",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/francoislangot/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([7.51739835, 3.00535761, 1.46651755])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_c   = (psi[1,0], psi[1,1], psi[1,2], 1, yy_c[1], sigma)\n",
    "pd_c0     = [6 , 2 , 2 ] \n",
    "calibsol  = fsolve(hetero_calib, pd_c0, args=param_c) \n",
    "calibsol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8945fb87",
   "metadata": {},
   "source": [
    "We check the solution: to find the health gradient with the model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "214a0e9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.        , 0.081768  , 0.18348814, 0.33634386]), 1.2759999999998748)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv = np.zeros(yy.size)\n",
    "\n",
    "for ii in range(yy.size):\n",
    "    param_m   = (1, yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolv[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "msolv, ( 1-np.exp(-(calibsol[2]+calibsol[0]*msolv[-1])) )/( 1-np.exp(-(calibsol[2]+calibsol[0]*msolv[0])) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1543e6dd",
   "metadata": {},
   "source": [
    "## Solution for the Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "986c98f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_s           = (psi[0,0], psi[0,1], yy_c[0], sigma, calibsol[1], calibsol[0])\n",
    "pd_c0             = [0.9 , calibsol[2]]\n",
    "calibsol_eu       = fsolve(hetero_calib_eu, pd_c0, args=param_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7a0199ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.56730391, 1.89073796])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calibsol_eu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d3a89f8",
   "metadata": {},
   "source": [
    "## Individual decisions and health gradient: Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "dd23a4f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.01549786, 0.09418598, 0.16101862, 0.25451866]), 1.1294746245585219)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv_eu = np.zeros(yy.size)\n",
    "for jj in range(yy.size):\n",
    "    param_m         = (calibsol_eu[0], yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eu[jj]    = max(0,fsolve(medexp, m0, args=param_m))\n",
    "grads_eu = ( 1-np.exp(-(calibsol_eu[1]+calibsol[0]*msolv_eu[-1])) )/(1-np.exp(-(calibsol_eu[1]+calibsol[0]*msolv_eu[0])))\n",
    "    \n",
    "msolv_eu, grads_eu\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "597eefdf",
   "metadata": {},
   "source": [
    "## Demand : price elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f3322eb8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.7497906679729238, 0.6677553518387034)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolvp = np.zeros(yy.size)\n",
    "msolvm = np.zeros(yy.size)\n",
    "\n",
    "for ii in range(yy.size):\n",
    "    param_m   = (1.1, yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvp[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m   = (0.9, yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvm[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "dpop    = (1.1-1)/1\n",
    "level0  = s_c[1]\n",
    "level1p = np.mean(msolvp)\n",
    "level1m = np.mean(msolvm)\n",
    "\n",
    "((level1m-level0)/level0)/dpop, -((level1p-level0)/level0 )/dpop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6a3d7e61",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.7087730099058136"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-(((level1m-level0)/level0)/dpop + (-((level1p-level0)/level0 )/dpop))/2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c4bbc61",
   "metadata": {},
   "source": [
    "## Demand : income elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "92110032",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.1746605940831996, 1.2636166918805645)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolvp = np.zeros(yy.size)\n",
    "msolvm = np.zeros(yy.size)\n",
    "\n",
    "for ii in range(yy.size):\n",
    "    param_m   = (1, 1.1*yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvp[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m   = (1, 0.9*yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvm[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "dyoy    = (1.1-1)/1\n",
    "level0  = s_c[1]\n",
    "level1p = np.mean(msolvp)\n",
    "level1m = np.mean(msolvm)\n",
    "\n",
    "((level1p-level0)/level0)/dyoy, -((level1m-level0)/level0)/dyoy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f964f2f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.219138642981882"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((level1p-level0)/level0)/dyoy)+(-((level1m-level0)/level0)/dyoy))/2 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3af07c7a",
   "metadata": {},
   "source": [
    "## Demand : price elasticity in EU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ba949ad1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.040449852797365, 0.9226845436459346)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv_eup = np.zeros(yy.size)\n",
    "msolv_eum = np.zeros(yy.size)\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m         = (calibsol_eu[0]*1.1, yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eup[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m         = (calibsol_eu[0]*0.9, yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eum[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "level0  = s_c[0]/(calibsol_eu[0]/y_c[0]) \n",
    "level1p = np.mean(msolv_eup) \n",
    "level1m = np.mean(msolv_eum) \n",
    "\n",
    "((level1m-level0)/level0)/dpop, -((level1p-level0)/level0)/dpop "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8e41e4a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.9815671982216497"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-((((level1m-level0)/level0)/dpop)+(-((level1p-level0)/level0)/dpop))/2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdaae432",
   "metadata": {},
   "source": [
    "## Demand : income elasticity in EU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "ee03d8bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.7356186895066263, 1.7690924707129616)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv_eup = np.zeros(yy.size)\n",
    "msolv_eum = np.zeros(yy.size)\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m         = (calibsol_eu[0], 1.1*yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eup[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m         = (calibsol_eu[0], 0.9*yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eum[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "level0  = s_c[0]/(calibsol_eu[0]/y_c[0])    \n",
    "level1p = np.mean(msolv_eup) \n",
    "level1m = np.mean(msolv_eum)\n",
    "\n",
    "((level1p-level0)/level0)/dyoy, -((level1m-level0)/level0)/dyoy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1e6bf63d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.752355580109794"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((level1p-level0)/level0)/dpop)+(-((level1m-level0)/level0)/dpop))/2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf13ddbd",
   "metadata": {},
   "source": [
    "# Forecasting of the GDP share of health expenditures\n",
    "\n",
    "Observation: the shift in GDP per capita\n",
    "\n",
    "Implication: the implicit shift in healthcare prices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5fa8267",
   "metadata": {},
   "source": [
    "### GDP per capita\n",
    "\n",
    "GDP share of health expenditures: OECD data\n",
    "\n",
    "Real GDP per capita: WB data (World Development Indicators)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "2d32641c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1970\n",
    "GDPpc_eu70 = [17887.10176179620, 26927.95903692180, 17519.42717022000, 15729.85644157820, 21354.25810368360, 24388.29821013090, 11431.59455874400]\n",
    "PMoY_eu70  = [5.713, 7.676666666666667,   5.184, float(\"NAN\"), 5.434246479956248,  5.506, 3.145]\n",
    "GDPpc70    = [np.nanmean(GDPpc_eu70), 25279.18879113590]\n",
    "PMoY70     = [np.nanmean(PMoY_eu70) , 6.232]\n",
    "# 1990\n",
    "GDPpc_eu90 = [29473.85438511120,   39295.24344249190,   28512.10870941340,   27479.57862963770,   31093.54486501490,   34570.30573399390,   18962.44605133250]\n",
    "PMoY_eu90  = [8.031, 8.038, 7.958, 7.015, 7.082, 7.256, 6.101]\n",
    "GDPpc90    = [np.nanmean(GDPpc_eu90) , 39278.59534978750]\n",
    "PMoY90     = [np.nanmean(PMoY_eu90) , 11.273]\n",
    "# 2000\n",
    "GDPpc_eu00 = [34476.20802831070, 49241.92800413450, 33583.85744208620, 32337.89674255510, 40440.67534426340, 41117.15243035039, 23928.67165633810]\n",
    "PMoY_eu00  = [9.827999999999999, 8.103999999999999, 9.541, 7.58, 7.056, 7.412, 6.8160]\n",
    "GDPpc00 = [np.nanmean(GDPpc_eu00) , 48689.03128985530]\n",
    "PMoY00  = [np.nanmean(PMoY_eu00) , 12.502]\n",
    "# 2017\n",
    "GDPpc_eu10 = [42622.40993268510, 55735.76490114850, 37678.92729753470, 31231.66444646460, 46978.44880120770, 52576.81028210510, 27213.45166906160]\n",
    "PMoY_eu10  = [11.272, 10.218, 11.458, 8.9010, 10.143, 10.916, 8.840999999999999]\n",
    "GDPpc10    = [np.nanmean(GDPpc_eu10) , 58387.77580835719]\n",
    "PMoY10     = [np.nanmean(PMoY_eu10) , 17.149999999999999]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "4e128ceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "psi0    = [i/100 for i in PMoY90]\n",
    "psi1    = [i/100 for i in PMoY10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba04b0a3",
   "metadata": {},
   "source": [
    "## Income heterogeneity by country (data)\n",
    "We rescale the income level and compute the quartiles in 1970 and 2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "aab01b00",
   "metadata": {},
   "outputs": [],
   "source": [
    "y70_c = y_c*np.array(GDPpc70)/np.array(GDPpc00)\n",
    "y90_c = y_c*np.array(GDPpc90)/np.array(GDPpc00)\n",
    "y00_c = y_c*np.array(GDPpc00)/np.array(GDPpc00)\n",
    "y10_c = y_c*np.array(GDPpc10)/np.array(GDPpc00)\n",
    "\n",
    "yy70_c = np.zeros((Ncount,yy.size))\n",
    "yy90_c = np.zeros((Ncount,yy.size))\n",
    "yy00_c = np.zeros((Ncount,yy.size))\n",
    "yy10_c = np.zeros((Ncount,yy.size))\n",
    "\n",
    "for ii in range(Ncount):\n",
    "    yy    = [norm(0,np.sqrt(sig_c[ii])).ppf(q) for q in qs]\n",
    "    yy    = np.exp(yy)\n",
    "    yy70_c[ii,:] = y70_c[ii]*yy/np.mean(yy)\n",
    "    yy90_c[ii,:] = y90_c[ii]*yy/np.mean(yy)\n",
    "    yy00_c[ii,:] = y00_c[ii]*yy/np.mean(yy)\n",
    "    yy10_c[ii,:] = y10_c[ii]*yy/np.mean(yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "de283a97",
   "metadata": {},
   "outputs": [],
   "source": [
    "prevsol0 = [0,0]\n",
    "prevsol1 = [0,0]\n",
    "#                  pm/y     income    sigma  phi          alpha_1         alpha_0\n",
    "param_s        = (psi0[0], yy90_c[0], sigma, calibsol[1], calibsol_eu[1], calibsol[0])\n",
    "pd_c0          = .1 #.3\n",
    "prevsol0[0]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n",
    "param_s        = (psi0[1], yy90_c[1], sigma, calibsol[1], calibsol[2], calibsol[0])\n",
    "pd_c0          = .8\n",
    "prevsol0[1]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n",
    "\n",
    "param_s        = (psi1[0], yy10_c[0], sigma, calibsol[1], calibsol_eu[1], calibsol[0])\n",
    "pd_c0          = .5#.5\n",
    "prevsol1[0]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n",
    "param_s        = (psi1[1], yy10_c[1], sigma, calibsol[1], calibsol[2], calibsol[0])\n",
    "pd_c0          = 1.2\n",
    "prevsol1[1]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "43bea0f2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([[0.23604486640425093], [0.4808111082647214]],\n",
       " [[0.5199008531699448], [1.3770087808513378]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prevsol0, prevsol1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "825acd69",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAATLElEQVR4nO3dYYhdZ37f8e9vZYtdSFvZ9dgISVRmGULVQLRCaAV5s8TZImlLxwvd1iqxhTHIohZsIKVMkxdx8srdZuNiEBJyV0SmS4QhKR6caYRR1oSFytF469Va6whPhbue9WBPHNaJcYiR/e+LOWpu717pnpHGmrWe7wcO95zn+T/3Pg9c7k/n6Ny5qSokSe35zFpPQJK0NgwASWqUASBJjTIAJKlRBoAkNeq2tZ7AStx11121devWtZ6GJH2qvPzyy39ZVRPD7Z+qANi6dStzc3NrPQ1J+lRJ8n9GtXsJSJIaZQBIUqMMAElqlAEgSY3qFQBJ9iS5mGQ+yfSI/iR5qus/n2THUP+6JP8ryfMDbXcmeSHJ693jHTe+HElSX2MDIMk64AiwF9gG7E+ybahsLzDZbQeBo0P9XwdeG2qbBs5U1SRwpjuWJN0kfc4AdgHzVXWpqj4ETgFTQzVTwDO17CywIclGgCSbga8A/3XEmJPd/kng/utbgiTpevQJgE3AmwPHC11b35r/AvwH4OOhMfdU1SJA93j3qBdPcjDJXJK5paWlHtOVJPXRJwAyom34RwRG1iT5F8A7VfXyimd25UmqjlfVzqraOTHxU19kkyRdpz7fBF4Atgwcbwbe6lnzr4B/mWQf8FngHyb5b1X1q8DbSTZW1WJ3ueid612EdCvYOv3Haz0F/Qx744mvrPpz9jkDOAdMJrk3yXrgAWBmqGYGeKi7G2g38F5VLVbVf6yqzVW1tRv3p92H/5UxB7r9A8BzN7oYSVJ/Y88AqupyksPAaWAdcKKqLiQ51PUfA2aBfcA88AHwcI/XfgJ4NskjwI+Ar13fEiRJ16PXH4OrqlmWP+QH244N7Bfw2JjneBF4ceD4XeC+/lOVJK0mvwksSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjeoVAEn2JLmYZD7J9Ij+JHmq6z+fZEfX/tkkf57k+0kuJPntgTGPJ/lxkle6bd/qLUuSNM7Yn4RMsg44AnwZWADOJZmpqh8OlO0FJrvti8DR7vHvgF+uqveT3A58N8n/qKqz3bgnq+p3V285kqS++pwB7ALmq+pSVX0InAKmhmqmgGdq2VlgQ5KN3fH7Xc3t3VarNXlJ0vXrEwCbgDcHjhe6tl41SdYleQV4B3ihql4aqDvcXTI6keSOUS+e5GCSuSRzS0tLPaYrSeqjTwBkRNvwv+KvWlNVH1XVdmAzsCvJL3T9R4HPA9uBReCbo168qo5X1c6q2jkxMdFjupKkPvoEwAKwZeB4M/DWSmuq6ifAi8Ce7vjtLhw+Bp5m+VKTJOkm6RMA54DJJPcmWQ88AMwM1cwAD3V3A+0G3quqxSQTSTYAJPkc8CvAX3THGwfGfxV49caWIklaibF3AVXV5SSHgdPAOuBEVV1IcqjrPwbMAvuAeeAD4OFu+EbgZHcn0WeAZ6vq+a7vG0m2s3yp6A3g0dValCRpvLEBAFBVsyx/yA+2HRvYL+CxEePOA1+4ynM+uKKZSpJWld8ElqRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEb1CoAke5JcTDKfZHpEf5I81fWfT7Kja/9skj9P8v0kF5L89sCYO5O8kOT17vGO1VuWJGmcsQHQ/Z7vEWAvsA3Yn2TbUNleYLLbDgJHu/a/A365qn4R2A7s6X40HmAaOFNVk8CZ7liSdJP0OQPYBcxX1aWq+hA4BUwN1UwBz9Sys8CGJBu74/e7mtu7rQbGnOz2TwL338A6JEkr1CcANgFvDhwvdG29apKsS/IK8A7wQlW91NXcU1WLAN3j3aNePMnBJHNJ5paWlnpMV5LUR58AyIi26ltTVR9V1XZgM7AryS+sZIJVdbyqdlbVzomJiZUMlSRdw209ahaALQPHm4G3VlpTVT9J8iKwB3gVeLu7TLSYZCPLZwifmK3Tf/xJPr0+xd544itrPQVpTfQ5AzgHTCa5N8l64AFgZqhmBniouxtoN/Be98E+kWQDQJLPAb8C/MXAmAPd/gHguRtbiiRpJcaeAVTV5SSHgdPAOuBEVV1IcqjrPwbMAvuAeeAD4OFu+EbgZHcn0WeAZ6vq+a7vCeDZJI8APwK+tnrLkiSN0+cSEFU1y/KH/GDbsYH9Ah4bMe488IWrPOe7wH0rmawkafX4TWBJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqVK8ASLInycUk80mmR/QnyVNd//kkO7r2LUm+k+S1JBeSfH1gzONJfpzklW7bt3rLkiSNM/YnIbvf8z0CfBlYAM4lmamqHw6U7QUmu+2LwNHu8TLw61X1vST/AHg5yQsDY5+sqt9dveVIkvrqcwawC5ivqktV9SFwCpgaqpkCnqllZ4ENSTZW1WJVfQ+gqv4GeA3YtIrzlyRdpz4BsAl4c+B4gZ/+EB9bk2Qryz8Q/9JA8+HuktGJJHf0nbQk6cb1CYCMaKuV1CT5OeAPgV+rqr/umo8Cnwe2A4vAN0e+eHIwyVySuaWlpR7TlST10ScAFoAtA8ebgbf61iS5neUP/29X1R9dKaiqt6vqo6r6GHia5UtNP6WqjlfVzqraOTEx0WO6kqQ++gTAOWAyyb1J1gMPADNDNTPAQ93dQLuB96pqMUmAbwGvVdXvDQ5IsnHg8KvAq9e9CknSio29C6iqLic5DJwG1gEnqupCkkNd/zFgFtgHzAMfAA93w38JeBD4QZJXurbfqKpZ4BtJtrN8qegN4NFVWpMkqYexAQDQfWDPDrUdG9gv4LER477L6P8foKoeXNFMJUmrym8CS1KjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqN6BUCSPUkuJplPMj2iP0me6vrPJ9nRtW9J8p0kryW5kOTrA2PuTPJCkte7xztWb1mSpHHGBkCSdcARYC+wDdifZNtQ2V5gstsOAke79svAr1fVPwV2A48NjJ0GzlTVJHCmO5Yk3SR9zgB2AfNVdamqPgROAVNDNVPAM7XsLLAhycaqWqyq7wFU1d8ArwGbBsac7PZPAvff2FIkSSvRJwA2AW8OHC/w9x/ivWuSbAW+ALzUNd1TVYsA3ePdvWctSbphfQIgI9pqJTVJfg74Q+DXquqv+08PkhxMMpdkbmlpaSVDJUnX0CcAFoAtA8ebgbf61iS5neUP/29X1R8N1LydZGNXsxF4Z9SLV9XxqtpZVTsnJiZ6TFeS1EefADgHTCa5N8l64AFgZqhmBniouxtoN/BeVS0mCfAt4LWq+r0RYw50+weA5657FZKkFbttXEFVXU5yGDgNrANOVNWFJIe6/mPALLAPmAc+AB7uhv8S8CDwgySvdG2/UVWzwBPAs0keAX4EfG3VViVJGmtsAAB0H9izQ23HBvYLeGzEuO8y+v8HqKp3gftWMllJ0urxm8CS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUqF4BkGRPkotJ5pNMj+hPkqe6/vNJdgz0nUjyTpJXh8Y8nuTHSV7ptn03vhxJUl9jAyDJOuAIsBfYBuxPsm2obC8w2W0HgaMDfb8P7LnK0z9ZVdu7bfYqNZKkT0CfM4BdwHxVXaqqD4FTwNRQzRTwTC07C2xIshGgqv4M+KvVnLQk6cb1CYBNwJsDxwtd20prRjncXTI6keSOHvWSpFXSJwAyoq2uo2bYUeDzwHZgEfjmyBdPDiaZSzK3tLQ05iklSX31CYAFYMvA8Wbgreuo+f9U1dtV9VFVfQw8zfKlplF1x6tqZ1XtnJiY6DFdSVIffQLgHDCZ5N4k64EHgJmhmhngoe5uoN3Ae1W1eK0nvfJ/BJ2vAq9erVaStPpuG1dQVZeTHAZOA+uAE1V1Icmhrv8YMAvsA+aBD4CHr4xP8gfAl4C7kiwAv1VV3wK+kWQ7y5eK3gAeXb1lSZLGGRsAAN0tmrNDbccG9gt47Cpj91+l/cH+05QkrTa/CSxJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVG9AiDJniQXk8wnmR7RnyRPdf3nk+wY6DuR5J0krw6NuTPJC0le7x7vuPHlSJL6GhsASdYBR4C9wDZgf5JtQ2V7gcluOwgcHej7fWDPiKeeBs5U1SRwpjuWJN0kfc4AdgHzVXWpqj4ETgFTQzVTwDO17CywIclGgKr6M+CvRjzvFHCy2z8J3H8d85ckXac+AbAJeHPgeKFrW2nNsHuqahGge7x7VFGSg0nmkswtLS31mK4kqY8+AZARbXUdNdelqo5X1c6q2jkxMbEaTylJol8ALABbBo43A29dR82wt69cJuoe3+kxF0nSKukTAOeAyST3JlkPPADMDNXMAA91dwPtBt67cnnnGmaAA93+AeC5FcxbknSDxgZAVV0GDgOngdeAZ6vqQpJDSQ51ZbPAJWAeeBr4d1fGJ/kD4H8CP59kIckjXdcTwJeTvA58uTuWJN0kt/UpqqpZlj/kB9uODewX8NhVxu6/Svu7wH29ZypJWlV+E1iSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIa1SsAkuxJcjHJfJLpEf1J8lTXfz7JjnFjkzye5MdJXum2fauzJElSH2MDIMk64AiwF9gG7E+ybahsLzDZbQeBoz3HPllV27ttFknSTdPnDGAXMF9Vl6rqQ+AUMDVUMwU8U8vOAhuSbOw5VpK0BvoEwCbgzYHjha6tT824sYe7S0Ynktwx6sWTHEwyl2RuaWmpx3QlSX30CYCMaKueNdcaexT4PLAdWAS+OerFq+p4Ve2sqp0TExM9pitJ6uO2HjULwJaB483AWz1r1l9tbFW9faUxydPA871nLUm6YX3OAM4Bk0nuTbIeeACYGaqZAR7q7gbaDbxXVYvXGtv9H8EVXwVevcG1SJJWYOwZQFVdTnIYOA2sA05U1YUkh7r+Y8AssA+YBz4AHr7W2O6pv5FkO8uXhN4AHl3FdUmSxuhzCYjuFs3ZobZjA/sFPNZ3bNf+4IpmKklaVX4TWJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhrVKwCS7ElyMcl8kukR/UnyVNd/PsmOcWOT3JnkhSSvd493rM6SJEl9jA2AJOuAI8BeYBuwP8m2obK9wGS3HQSO9hg7DZypqkngTHcsSbpJ+pwB7ALmq+pSVX0InAKmhmqmgGdq2VlgQ5KNY8ZOASe7/ZPA/Te2FEnSSvT5UfhNwJsDxwvAF3vUbBoz9p6qWgSoqsUkd4968SQHWT6rAHg/ycUec9Z4dwF/udaT+FmQ/7TWM9BV+B4dcIPv038yqrFPAGREW/Ws6TP2mqrqOHB8JWM0XpK5qtq51vOQrsb36CevzyWgBWDLwPFm4K2eNdca+3Z3mYju8Z3+05Yk3ag+AXAOmExyb5L1wAPAzFDNDPBQdzfQbuC97vLOtcbOAAe6/QPAcze4FknSCoy9BFRVl5McBk4D64ATVXUhyaGu/xgwC+wD5oEPgIevNbZ76ieAZ5M8AvwI+NqqrkzjeFlNP+t8j37CUrWiS/KSpFuE3wSWpEYZAJLUKAPgFpfkoySvDGzTXfsbSe4aqPtSkufXbqZqVZKtSV4dans8yb9PsjvJS91797Ukj6/RNG9Jfb4HoE+3v62q7Ws9Cek6nQT+dVV9v/vTMj+/1hO6lRgAkn6W3Q1c+YsBHwE/XNvp3Fq8BHTr+9zQJaB/s9YTklbgSeBikv+e5NEkn13rCd1KPAO49V3tEtCo+3+9J1hr4Wrvu6qq30nybeCfA/8W2A986WZN7FbnGUC73gUGf4PhTvzDW1obw+9FGHg/VtX/rqqjwH3ALyb5xzd5frcsA6BdLwIPwv/73YZfBb6zlhNSm6rqfWAxyX2w/GNRwB7gu0m+kuTKH5WcBD4CfrImE70F+U3gW1ySj4AfDDT9SVVNJ/lHLP9wzz9j+a+2/gkwXVUfr8E01bjuh6KO8PdnAv+5qr6d5BSwg+U/MXMZ+M2qOr1G07zlGACS1CgvAUlSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1Kj/C5c9BVkZYb4KAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "columns = ('EU', 'US')\n",
    "Price0 = sum(prevsol0,[])\n",
    "Price1 = sum(prevsol1,[]) \n",
    "gPrice = [((Price1[ii]/Price0[ii])**(1/(2017-1990)))-1 for ii in np.arange(len(columns))]\n",
    "plt.bar(np.arange(len(columns)), gPrice)\n",
    "plt.xticks(np.arange(len(columns)), labels=columns)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9bd05d30",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.029676881213976047, 0.039739460012813366]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gPrice"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23335f30",
   "metadata": {},
   "source": [
    "### Size of the error when price changes are omitted\n",
    "\n",
    "If prices change, then the model predicts exactly the GDP share of health expenditures (targeted data)\n",
    "\n",
    "$\\Rightarrow$ If prices do not change, then it exists a gap between model and data: this gap gives the error if price change are neglected "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "a119d338",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([8.112481459718433, 13.91943187305987],\n",
       " [7.3544285714285715, 11.273],\n",
       " [10.44032415091521, 15.510530852904944],\n",
       " [10.249857142857142, 17.15])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv0_90 = np.zeros((2,yy.size))\n",
    "pmoys0_90 = [0,0]\n",
    "msolv0_10 = np.zeros((2,yy.size))\n",
    "pmoys0_10 = [0,0]\n",
    "\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m        = (calibsol_eu[0], yy90_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0             = .1\n",
    "    msolv0_90[0,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m        = (calibsol_eu[0], yy10_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0             = .1\n",
    "    msolv0_10[0,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m        = (1, yy90_c[1,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0             = .1\n",
    "    msolv0_90[1,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m        = (1, yy10_c[1,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0             = .1\n",
    "    msolv0_10[1,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "price_bench = [calibsol_eu[0],1] \n",
    "for ii in range(2):            \n",
    "    pmoys0_90[ii]= 100*price_bench[ii]*np.mean(msolv0_90[ii,:])/np.mean(yy90_c[ii,:])\n",
    "    pmoys0_10[ii]= 100*price_bench[ii]*np.mean(msolv0_10[ii,:])/np.mean(yy10_c[ii,:])\n",
    "\n",
    "pmoys0_90, PMoY90, pmoys0_10, PMoY10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83f75020",
   "metadata": {},
   "source": [
    "Without price change, the model over-estimate the GDP share of health expenditures in the 1990s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "2cdc6ddb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.7580528882898614, 2.6464318730598695]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[pmoys0_90[ii]-PMoY90[ii] for ii in np.arange(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "267e1682",
   "metadata": {},
   "source": [
    "Without price changes, the model under-estimate the GDP share of health expenditures in the 2010s for the U.S."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "ef97b436",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.19046700805806793, -1.6394691470950544]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[pmoys0_10[ii]-PMoY10[ii] for ii in np.arange(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef478f4c",
   "metadata": {},
   "source": [
    "# Forcasting long run change in GDP share of Health Expenditures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7df75ef3",
   "metadata": {},
   "source": [
    "### Long run changes in GDP per capita and income dispersion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "cc424891",
   "metadata": {},
   "outputs": [],
   "source": [
    "gGDPpc    = [((GDPpc10[ii]/GDPpc90[ii])**(1/(2017-1990))-1) for ii in np.arange(2)]\n",
    "agGDPpc   = np.array(gGDPpc)\n",
    "agGDPpc10 = np.array(GDPpc10)\n",
    "\n",
    "horizon = 4\n",
    "fGDPpc  = np.zeros((2,horizon))\n",
    "fy_c    = np.zeros((2,horizon))\n",
    "fyy_c   = np.zeros((2,horizon,yy.size))\n",
    "\n",
    "for jj in range(2):\n",
    "    for ii in range(horizon):\n",
    "        fGDPpc[jj,ii] = (GDPpc10[jj]*(1+gGDPpc[jj])**(2030+(ii*10)-2017))/GDPpc00[jj]\n",
    "        fy_c[jj,ii] = y_c[jj]*fGDPpc[jj,ii]\n",
    "        yy    = [norm(0,np.sqrt(sig_c[jj])).ppf(q) for q in qs]\n",
    "        yy    = np.exp(yy)\n",
    "        fyy_c[jj,ii,:] = fy_c[jj,ii]*yy/np.mean(yy)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "2d40c0d9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.35719765, 1.53905893, 1.74528919, 1.97915382],\n",
       "       [1.45139469, 1.68093419, 1.94677557, 2.25466004]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fGDPpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "2b8441b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.01265430880932894, 0.014790792069827585]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gGDPpc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a307dd87",
   "metadata": {},
   "source": [
    "### Long run changes in prices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "0f77bc8f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.76038411, 1.01869143, 1.36474739, 1.82836075],\n",
       "       [2.28536491, 3.37443308, 4.98248599, 7.35684069]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fPrice  = np.zeros((2,horizon))\n",
    "\n",
    "for jj in range(2):\n",
    "    for ii in range(horizon):\n",
    "        fPrice[jj,ii] = (Price1[jj]*(1+gPrice[jj])**(2030+(ii*10)-2017))\n",
    "        \n",
    "fPrice"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e29d255",
   "metadata": {},
   "source": [
    "Europe with changes in price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "bbfbdab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx     = np.zeros((horizon,yy.size))\n",
    "yy_eu  = np.zeros(horizon)\n",
    "yy_eu0 = np.zeros(horizon)\n",
    "yy_us  = np.zeros(horizon)\n",
    "yy_us0 = np.zeros(horizon)\n",
    "\n",
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (fPrice[0,ii], fyy_c[0,ii,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "        m0               = .05#.01\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_eu[ii] = 100*fPrice[0,ii]*np.mean(xx[ii,:])/np.mean(fyy_c[0,ii,:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c852a7b4",
   "metadata": {},
   "source": [
    "Europe without price changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "ce3c0840",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (Price1[0], fyy_c[0,ii,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "        m0               = .01#.005\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_eu0[ii] = 100*Price1[0]*np.mean(xx[ii,:])/np.mean(fyy_c[0,ii,:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99c5497f",
   "metadata": {},
   "source": [
    "The U.S. with price changes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ce7da493",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (fPrice[1,ii], fyy_c[1,ii,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "        m0               = .01#.005\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_us[ii] = 100*fPrice[1,ii]*np.mean(xx[ii,:])/np.mean(fyy_c[1,ii,:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2079c7d",
   "metadata": {},
   "source": [
    "The U.S. without price changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "58fb5606",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = np.zeros((horizon,yy.size))\n",
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (Price1[1], fyy_c[1,ii,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "        m0               = .01#.005\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_us0[ii] = 100*Price1[1]*np.mean(xx[ii,:])/np.mean(fyy_c[1,ii,:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "8003a56b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6cklEQVR4nO3de3xU5bXw8d+TCyQk3IMICRBJSCSAgIKoQURBiKCk4inKsSpSizesbVELWou1r6dotX17SqsicryjvIdAgGIAES9BRC6iIiEQIEhIBBKuIQlJZtb7x0ymk2Qm15nMZLK+n08+M7Nn7z2LTbLXPPt59rOMiKCUUkrVFOTrAJRSSvknTRBKKaVc0gShlFLKJU0QSimlXNIEoZRSyqUQXwfgSVFRURIbG+vrMJRSqtXYsWNHoYj0cPVeQCWI2NhYtm/f7uswlFKq1TDGHHb3nl5iUkop5ZImCKWUUi5pglBKKeVSQPVBuFJRUUFeXh5lZWW+DiUghYWFERMTQ2hoqK9DUUp5WMAniLy8PDp27EhsbCzGGF+HE1BEhKKiIvLy8rjkkkt8HY5SysO8liCMMX2At4CLASuwSET+Zoz5AEi0r9YFOC0iw1xsnwucAyxApYiMaEocZWVlmhy8xBhD9+7dOXHihK9DUarNSU5OprCwsNbyqKgoNm/e7JHP8GYLohKYIyI7jTEdgR3GmA0icnvVCsaYl4AzdezjehGpfQQaSZOD9+ixVco3XCWHupY3hdcShIgUAAX25+eMMVlANLAHwNjOLNOAG7wVg1JKqaZrkVFMxphYYDiw1WnxtcAxEdnvZjMB1htjdhhjZnk5RKWUahXOnz/P2rVrW+SzvN5JbYyJBJYDvxKRs05vTQeW1rFpsojkG2MuAjYYY/aKyGcu9j8LmAXQt29fD0aulFL+obi4mE8++YSMjAw+++wzLly40CKf69UWhDEmFFtyeFdE0pyWhwBTgQ/cbSsi+fbH48AK4Eo36y0SkREiMqJHD5fTifiF4OBghg0b5vhZsGABubm5DB48uNp6zzzzDC+++GK9+1uxYgXGGPbu3eutkMnIyCAxMZH4+HgWLFjgtc9RStVWXFzMqlWreOihh7j66quZM2cO33zzDdOmTeOdd95pkRi8OYrJAK8DWSLylxpvjwf2ikiem20jgCB730UEMAF41luxOsvMzGTZsmUUFhYSFRXFtGnTGD16dLP3Gx4ezq5du6oty83NbfL+li5dyogRI3j//fd55plnmhWbxWIhPz+fPn36VFv28MMPs2HDBmJiYhg5ciRTpkwhKSmpWZ+llHLv7NmzfPzxx2RkZJCZmUlFRQU9e/bkjjvuICUlheHDhxMUZPteHxUV5XYUk6d48xJTMnAX8J0xZpd92ZMisha4gxqXl4wxvYHFIjIJ6AmssI+QCQHeE5EML8YK2JLD4sWLKS8vB2yjARYvXgzgkSThKcXFxXz66ads2LCBn/70p9USxB133IGIkJuby48//sg///lPJk+eXOf+5syZQ1xcHI888ohj2VdffUV8fDz9+/d37Dc9PV0ThFIedubMGUdS2Lx5MxUVFfTq1Ys777yTlJQUhg4d6kgKzjw1lLUu3hzFlAm4HAMpIjNcLMsHJtmfHwSGejqmt956i8OH3U5cyP79+6msrKy2rLy8nEWLFrFp0yaX2/Tr14+777673s8uLS1l2LBhjtfz5s1j1KhRDQu8hpUrVzJ+/Hguu+wyIiIi2LlzJ5dffjkA33zzDT/5yU/44IMPyMzM5De/+U2tBDF+/Hh+/PFHwHazW1ZWFklJSfTt25fU1FQAjh49Wq1FERMTw9atW1FKNd/p06fZuHEjGRkZbNmyhYqKCqKjo7nrrrtISUlhyJAhLpNCSwv4O6kbo2ZyqG95Y7i6xOQuWdV3b8HSpUuZNcs2sGvatGksXbqUyy+/nNLSUgoLC5k/fz4ASUlJnDp1qtb2H330EQBffPEFTzzxBCUlJYSFhVVbR0QaHZdSyr1Tp07x0UcfkZGRwZdffkllZSXR0dHcfffdjqTgb39jbSpB1PdN/5e//KXba3pPP/20x+Pp3r17rRP4yZMn65y2oqioiK+++oq0NFuf/+233851113HCy+8wO7duxkwYIDjZL9z506GDq3dEKtqQRw5coROnToxYoTtJvXnnnvO0YKIiYnhyJEjjm3y8vLo3bt38/7BSrUxJ0+eZMOGDaxbt44vv/wSi8VCnz59uPfee0lJSWHQoEF+lxSctakEUZ9p06ZV64MAaNeuHdOmTfPK50VGRtKrVy82btzIuHHjOHnyJBkZGTz66KMAjBs3jrfeeovo6GjHNv/7v//LpEmTaN++PQCXXHIJF198MZmZmWRnZ/PDDz9QVlaGxWJh/vz5vPDCC7U+d/Xq1dx44428/vrrbi9zjRw5kv3793Po0CGio6N5//33ee+997xwFJQKLEVFRaxfv55169bx1VdfYbFY6NevH/fddx8pKSkMHDjQr5OCM00QTqo6opctW0ZRURHdu3f32Cimmn0QKSkpLFiwgLfeeouHH36YOXPmADB//nzi4uKwWq3k5OTQrVu3avtZunQp3377Lc6lVYuKinjvvfcICQnhzjvvZOzYsZw9e5Ynn3yS5ORkl/E899xzdfaBhISEsHDhQiZOnIjFYmHmzJkMGjSo6QdAqQB24sQJR1LYtm0bVquV2NhYfvGLX3DTTTeRmJjYapKCM+PqWnNrNWLECKlZcjQrK4uBAwf6KKKm2717N0uWLOEvf6k5Qti9MWPG8Nprr5GYmFj/yh7UWo+xUs1x/PjxaklBROjfvz8pKSmkpKSQkJDQKpKCMWaHu8lQtQXhpwYPHtyo5ABw4MABBgwY4KWIlFLHjh1j3bp1rFu3jh07diAixMfH8/DDD5OSkhJwf3+aIALI0aNHfR2CUgGnoKCAdevWkZGRwddffw1AQkICs2fPJiUlhfj4eB9H6D2aIJRSqob8/HxHUqganp6YmMijjz7KxIkTiYuL822ALUQThFJKYRvKXZUUvv32WwAGDhzIr3/9ayZOnNgmqyZqglBKtVlHjhwhIyODjIwMdu/eDcCgQYOYM2cOEydOpF+/fj6O0Lc0QSil2pTDhw87Wgrff/89AEOGDOGxxx4jJSWl2hQzbZ0mCKVUwDt06JAjKWRlZQEwdOhQnnjiCSZOnEhMTIyPI/RPmiCUUgHpwIEDjqSQnZ0NwPDhw5k7dy4TJ07UqWMaQBOEUipg5OTkkJGRwbp169i3bx8Al19+OfPmzWPixIn06tXLxxG2LpoglFKtloiwf/9+R0shJycHYwxXXHEFTz31FBMnTqRnz56+DrPV8v2E421AfaVFn3vuOQYNGsRll13GsGHDGlx3wdtlR7XkqPJHIsLevXv529/+xqRJk7jlllv4xz/+Qbdu3Xj66af59NNPeffdd7n77rs1OTSTtiCcJCcnu53u21vVm7Zs2cKaNWvYuXMn7du3p7CwsNpssnXxVNlRLTmq/EFdf3+ZmZns3bvXMSQ1NzeXoKAgRo4cyc9+9jMmTJiAP9ekb600QThx9ctZ13JPKCgoICoqyjF9d0PryborO6olR1VrVdff38SJEzl8+DBBQUGMGjWKGTNmcOONN3q0/rKqrU0liOeee67Jl2Puuusul8svvfRSnnrqqSbHNGHCBJ599lkSEhIYP368owBQfdyVHdWSoyoQxcTEcN999zF+/PhaU+Ar72lTCcJX3E35a4whMjKSHTt28Pnnn7Np0yZuv/12FixYwIwZM+rcp6uyowMHDtSSo6pVqvrC4s6SJUtaKBLlzGsJwhjTB3gLuBiwAotE5G/GmGeAXwAn7Ks+KSJrXWyfAvwNCAYWi0ize0nr+6ZfVx2Ft99+u8mfW19p0eDgYMaOHcvYsWMZMmQIb775Zp0Jwl3Z0WnTpmnJUdVqFBcXs2HDBtLT0/nyyy99HY5ywZujmCqBOSIyELgKeNgYU3UB+68iMsz+4yo5BAP/AG4CkoDpTtu2Os6lRQFHadHRo0eTnZ3N/v37Hevu2rXLMf/LuHHjXE7h7a7saEZGhqPk6Pnz55k/fz6//vWva22/evVqunTpwvr16zly5Ai7d+9m9+7djuQA1UuOlpeX8/777zNlyhSPHhfV9lRWVvL555/z2GOPMXr0aObOnUteXh4PP/ywr0NTLnitBSEiBUCB/fk5Y0wWEF33Vg5XAjkichDAGPM+kArs8UasVaKiotyOomgud6VFd+zYwSOPPMLp06cJCQkhPj6eRYsWuS05Cu7Ljl5xxRVaclT5pb1797Jy5UrWrFnDiRMn6NSpE1OmTCE1NZXLL78cYwzvv/++1/7+VNO0SMlRY0ws8BkwGPgNMAM4C2zH1so4VWP9/wBSROQ+++u7gFEiMtvFvmcBswD69u17xeHDh6u931rLYWrJUdXaHTt2jNWrV5Oens6+ffsIDQ1lzJgxpKamcv3119OuXTtfh6jwcclRY0wksBz4lYicNca8DPwREPvjS8DMmpu52JXLTCYii4BFYKtJ7am4fU1LjqrW6Pz582zYsIFVq1bxxRdfICIMGzaM3//+90yaNImuXbv6OkTVCF5NEMaYUGzJ4V0RSQMQkWNO778GrHGxaR7gPOduDJDvxVADgpYcVb5gsVj48ssvWblyJR999BElJSXExMTw4IMPMmXKlDZZaCdQeHMUkwFeB7JE5C9Oy3vZ+ycAbgV2u9h8GzDAGHMJcBS4A/hPb8WqlGq8vXv3kp6ezpo1azh+/DgdO3bk5ptvJjU1lSuuuEKHRQcAb7YgkoG7gO+MMbvsy57ENiJpGLZLRrnA/QDGmN7YhrNOEpFKY8xsYB22Ya5LROR7L8aqlGqA48ePs2bNGtLT09m7dy8hISHV+hWqRtapwODNUUyZuO5LqDWs1b5+PjDJ6fVad+sqpVpOSUkJH330Eenp6XzxxRdYrVYuu+wynn76aSZNmqR3NgcwvZNaKVWLxWJh69atpKens379ekpKSoiOjmbWrFmkpqY65uhSgU0ThFLKYd++faSnp7N69WqOHTtGZGQkkydPdvQrBAVphYC2RBOEUm3ciRMn+Ne//sXKlSvJysoiJCSEa6+9lrlz53LDDTfUmqdLtR2aIJRqg0pLSx39Cps3b8ZqtTJ48GB+97vfMXnyZO1XUIAmCKXaDKvVytatW1m1ahUZGRmUlJTQq1cvZs2axZQpU4iLi/N1iMrPaIJoIcHBwQwZMsTx+o477uCOO+7g5ptvZvfuf98K8swzzxAZGcljjz1W5/5WrFjB1KlTycrK4tJLL/VKzBkZGTz66KNYLBbuu+8+5s6d65XPUd6Vk5Pj6FcoKCggIiKCm266idTUVEaOHKn9CsotTRAuXLhwgT179pCUlOSxcd3h4eHs2rWr2rLc3Nwm789T5Ubd0ZKjrVthYSFr1qxh1apVfP/99wQHBzN69Ggef/xxbrjhBsLDw30domoFNEG4kJuby5kzZ8jNzW3xie8awl25UWhayVFXtORo61NWVsbGjRtJT08nMzMTi8XCoEGDePLJJ5k8ebLOiqoarU0liP3791NcXOz2/TNnzlR7XVBQQEGBbVaQzp07u9wmMjKyQRPklZaWMmzYMMfrefPm1Tnddl3clRsFGlRy9Nprr+XcuXO19vviiy8yfvx4QEuOthZWq5Vt27aRnp5ORkYG58+f5+KLL+bnP/85qampxMfH+zpE1Yq1qQRRn44dO1JWVkZFRYVjWWhoqEeG+bm6xFRzavIq9c1h46rc6OWXX05paWmDSo5+/vnn9carJUf924EDB1i1ahWrVq0iPz+fDh06kJKSQmpqKldeeaX2KyiPaFMJoiHf9LOzsykoKCAoKAir1UpUVJTXLjPVV4rUFXflRl944QV2797doJKjDWlBaMlR/1NUVMS//vUv0tPT2b17N0FBQSQnJzNnzhzGjRun/QrK49pUgmiIiooKevfuTe/evcnPz6e8vNxrn+VcinTcuHGOUqSPPvooYCs5+tZbbxEd/e9CfO7KjWZmZpKdne0oOWqxWJg/fz4vvPBCrc9tSAvCueRodHQ077//Pu+9956H/uWqocrKyvj4449ZtWoVn332GRaLhaSkJObNm8fkyZPp0aOHr0NUAUwTRA2DBw92PE9ISPDYfmv2QaSkpLBgwQK3pUjdlRx1V270vffeIyQkpMElR+ujJUd9x2q1smPHDlauXElGRgbFxcX07NmTe++9l9TUVI/+XipVF00QLcRisbhcnpSUxKZNm2ot37NnD7fddlutywaffPKJ28+oKjn6/PPPNyvWKpMmTWLSpEn1r6g84uDBg477FY4ePUqHDh2YOHGio18hODjY1yGqNkYThJ/SkqNtw8mTJ1m7di3p6el8++23BAUFcc011/CrX/2K8ePH06FDB1+HqNowTRABREuO+pfk5GQKCwtrLe/evTu///3vWblyJZ9//jmVlZVceuml/Pa3v+Xmm2/moosu8kG0StWmCUIpL3GVHMDWZ/Too4/So0cP7rnnHlJTU/3yhkylNEEo5QNLlizhqquu0n4F5dc0QSjlYSUlJaxfv77OdZo6ukyplqQJQikPEBF27txJWloaH374IefPn/d1SEo1m9cShDGmD/AWcDFgBRaJyN+MMX8GbgHKgQPAvSJy2sX2ucA5wAJUisgIb8WqVFMdO3aMlStXkpaWRm5urmPKi9tuu40777zT1+Ep1SzebEFUAnNEZKcxpiOwwxizAdgAzBORSmPM88A84Ldu9nG9iLju6VPKRy5cuMDGjRtJS0tzVGMbOXIk999/PxMnTiQiIgKAqKgolx3VOquqai28liBEpAAosD8/Z4zJAqJFxPni7JfAf3grBqU8RUT4/vvvSUtLY82aNZw5c4ZevXpx//33M3XqVPr27Vtrm82bN/sgUqU8p0X6IIwxscBwoOZ80TOBD9xsJsB6Y4wAr4rIIjf7ngXMAlz+kSrVHEVFRaxevZrly5ezb98+2rdvz4033sjUqVN1FJIKeF6fE9gYEwksB34lImedlj+F7TLUu242TRaRy4GbgIeNMWNcrSQii0RkhIiM8NTEZdnZ2cybN4/s7GyP7C83N7faHE9gKy364osvAvDcc88xaNAgLrvsMoYNG9bgugsrVqzAGMPevXs9EmdNGRkZJCYmEh8fz4IFC7zyGf6ooqKCjRs38vDDDzNmzBj+9Kc/ERYWxjPPPENmZiYvvfQSycnJmhxUwPNqC8IYE4otObwrImlOy+8BbgbGiavCA4CI5NsfjxtjVgBXAp95M16wJYeFCxdSXl7OwoULmT17tldvYtqyZQtr1qxh586dtG/fnsLCwgbPIOvNsqNtseTo/v37SUtLIz09naKiIqKiorj77ruZOnWqTmGi2iRvjmIywOtAloj8xWl5CrZO6etEpMTNthFAkL3vIgKYADzrrVirOCcHoEWSREFBAVFRUY7puxvagemu7KiWHG2cM2fO8K9//Yu0tDS+++47QkJCuP7665k6dSrXXnstoaGhvg5RKZ/xZgsiGbgL+M4Ys8u+7Engv4H2wAZ7hbIvReQBY0xvYLGITAJ6Aivs74cA74lIRnMD+uCDD8jLy3P5XklJCUePHq1VSa28vJy//vWvREdHu5w4LSYmhttvv73JMU2YMIFnn32WhIQExo8f7ygAVB93ZUe15Gj9LBYLW7ZsIS0tjQ0bNlBeXk5iYiLz5s1jypQptaZYV6qt8uYopkzAVY3KtW7Wzwcm2Z8fBGqXQvOiH3/80WWZTXs8/Pjjj45v043lrlSnMYbIyEh27NjB559/zqZNm7j99ttZsGABM2bMqHOfrsqODhw4UEuO1uHw4cOkpaWxcuVKfvzxR7p06cK0adOYOnUqSUlJrf7fp5Sntak7qev6pl/z8pKzdu3aNesyU32lRYODgxk7dixjx45lyJAhvPnmm3UmCHdlR6dNm6YlR2soLi4mIyODFStWsH37doKCghg9ejRz585l3LhxtGvXztchKuW32lSCqEtiYiKzZ8+ulSSamxyg7tKi2dnZBAUFOTpBd+3aRb9+/QDXJUfBfdnRjIwMLTmKrfWzfft20tLSyMjIoKSkhNjYWObMmUNqaio9e/b0dYhKtQqaIJzUTBKeSA5V3JUW3bFjB4888ginT58mJCSE+Ph4Fi1a5LbkKLgvO3rFFVe06ZKj+fn5rFy5khUrVvDDDz8QERHB5MmTmTp1KsOHD9dLSEo1knF33b01GjFihGzfvr3asqysLAYOHNio/WRnZ/PGG28wY8YMn83Tv3v3bpYsWdKoqnJVJUdbOuamHGNPKSsrY8OGDaSlpbFlyxZEhFGjRjF16lQmTJigFdmUqocxZoe7ue40QQSQ6Ohojhw5QlCQ1+9/rKalj7GI8O2337J8+XLWrl3LuXPniI6O5tZbb+UnP/lJtdFXSqm61ZUg9BJTAAn0kqMnTpwgPT2dFStWkJOTQ1hYGBMmTGDq1KmMGjWqxROjUoFOE4Tya+Xl5XzyySekpaXx2WefYbFYGD58OH/84x+56aab6Nixo69DVCpgaYJQfmnv3r0sX76c1atXc+rUKXr06MHMmTO59dZbiYuL83V4SrUJbSJBiIiOYPEST/ZhnTp1ijVr1pCWlsaePXsIDQ3lhhtu4LbbbiM5OZmQkDbx66qU3wj4v7iwsDCKioro3r27JgkPExGKioocN+Y1RWVlJZs3b2b58uV8/PHHVFRUkJSUxO9+9ztuvvlmunbt6sGIlVKNEfAJIiYmhry8PE6cOOHrUAJSWFgYMTExjd7u4MGDjplTjx8/TteuXZk+fTq33XYbl156qRciVUo1VsAniNDQUMeUFsq3iouLWbt2LWlpaXz99dcEBwczZswYnn76acaOHavTXijlZwI+QSjfslqtbN26lbS0NNavX09ZWRlxcXE8/vjjpKam4qkiT0opz9MEobwiLy+PFStWsGLFCo4ePUrHjh1JTU3ltttu47LLLtP+IKVaAU0QymNKS0tZt24daWlpbN26FWMMV199Nb/+9a+58cYbm9WZrZRqefUmCGPMCOBaoDdQCuwGPhKRk16OTfmZ5ORkCgsLay3v3LkzEyZMYO3atZw/f54+ffrwy1/+kltvvbVVThGulLJxmyCMMTOAXwKHgB1ANhAGjAZ+a4zZDTwtIj+0QJzKD7hKDmAr27lmzRpSUlKYOnUqI0aM0GkvlAoAdbUgIoBkESl19aYxZhgwANAEocjMzCQyMtLXYSilPMhtghCRf9S1oYjs8ng0ym999dVXdb6vyUGpwNPg6wDGmFuMMVuNMbuMMQ81YP0+xphNxpgsY8z3xphH7cu7GWM2GGP22x9d3iprjEkxxmQbY3KMMXMb/k9SniIibN68mTvvvJO77rrL1+EopVqY2wRhjKlZzPgu4CrgcuDBBuy7EpgjIgPt2z1sjEkC5gIbRWQAsNH+uuZnBwP/AG4CkoDp9m1VCxARNm3axO23387MmTM5cuQITz31lK/DUkq1sLr6IB4ytsHqvxeRH4EjwHOAFcivb8ciUgAU2J+fM8ZkAdFAKjDWvtqbwCfAb2tsfiWQIyIHAYwx79u329Ogf5VqEqvVyoYNG3j55ZfJysoiOjqaP/zhD0ydOpV27drx6quvuuyojoqK8kG0Silvq6sP4n57K+JVY8x24GngGqAD8MfGfIgxJhYYDmwFetqTByJSYIy5yMUm0dgSUpU8YJSbfc8CZgH07du3MWEpO4vFwtq1a3nllVfIyckhNjaWP/3pT9xyyy2EhoY61tu8ebMPo1RKtbQ674MQkW+AVGPMLcAq4E0RebsxH2CMiQSWA78SkbMNvIPW1Uou55UWkUXAIrCVHG1MbG1dRUUFq1atYtGiReTm5jJgwABeeuklbrrpJoKDg30dnlLKx+q6D+IB4H5sJ+YXgBRsl53WAf9HRD6vb+fGmFBsyeFdEUmzLz5mjOllbz30Ao672DQPcC4sHEMDLmuphikvL2f58uW89tprHD16lKSkJP7+978zfvx4vX9BKeVQZx+EiFxmjGkHbBGR94H/Nsa8je1yU50Jwt5/8TqQJSJ/cXprFXAPsMD+mO5i823AAGPMJcBR4A7gPxv4b1JulJaWsmzZMhYvXszx48cZOnQov//977nuuut0biSlVC11JYijxpg/AuHA3qqFInIK+E0D9p2MbeTTd8aYXfZlT2JLDMuMMT/HdpPdTwGMMb2BxSIySUQqjTGzgXVAMLBERL5v1L9MORQXF7N06VL+53/+h6KiIkaOHMnzzz/P1VdfrYlBKeWWcVcy0t5ymAhUABtExNKSgTXFiBEjZPv27b4Ow2+cPXuWd955hzfffJPTp0+TnJzMgw8+yMiRI30dmlLKTxhjdojICFfv1dWC6C0iq+vYqQGiRSSvuQEqzzp16hRvvvkmb7/9NsXFxVx//fU88MADDBs2zNehKaVakboSxJ+NMUHY+gh2ACewTdYXD1wPjAPmY+tQVn6gsLCQJUuWsHTpUkpKSpg4cSIPPPAASUl6j6FSqvHqug/ip/a7l+8EZgK9gBIgC1gLPCciZS0SparTjz/+yOLFi1m2bBkVFRVMmjSJBx54gAEDBvg6NKVUK1bffRB7AJ1jwU8dOXKE1157jbS0NESEKVOmMGvWLK3BrZTyCK0o1wrl5uby6quvkp6eTlBQEFOnTmXWrFnExMT4OjSlVADRBNGK7N+/n1deeYW1a9cSGhrKnXfeyc9//nMuvvhiX4emlApAmiBagT179vDyyy+zfv16OnTowL333svMmTN1kjyllFc1pCb1cmAJ8KGIWL0fkqryzTff8PLLL7Np0yYiIyN58MEHueeee+ja1WUJDaWU8qiGtCBeBu7FNs3G/wPeEJG99WyjmmHbtm28/PLLbN68mS5duvDoo4/ys5/9jE6dOvk6NKVUG1JvghCRj4CPjDGdgenABmPMEeA14B0RqfByjG2CiLBlyxb++c9/sm3bNrp3787jjz/O9OnTiYiI8HV4Sqk2qEF9EMaY7sDPsM2t9DXwLjAa22R7Y70VXFsgInz66af885//5JtvvqFnz5489dRT/PSnPyU8PNzX4SnV4rKzs3njjTeYMWMGiYmJvg6nTat3bmdjTBq2mVs7ALeIyBQR+UBEHgG0Un0TWa1W1q1bx6233sr9999PYWEhf/jDH/joo4+4++67NTkEkOzsbObNm0d2dravQ/F72dnZLFy4kJMnT7Jw4UI9Zj5W12R9VTUbbhCRj1s4riZpDZP1WSwWPvzwQ1555RX2799PbGws999/f63qbSowVJ3wysvLadeuHbNnz9ZvxW44H6sqesy8r67J+upKEB8CXbHVjM4AMkWk0ltBeoI/J4iKigpWr17Nq6++Sm5uLvHx8Tz44INavS2ANfWEJyKOH6vV6nh097yuZfU91resKftsyn5OnTrF3r17sVprD5QMCgpixIgRXHzxxYSGhrr8CQkJITQ0lHbt2jme13w/UIthNfeSXJMShH3DMGx9DDdhq+/wA7ZkkSEiPzQ6Ei/zxwRRXl5OWloaixYt4ujRowwcOJAHH3yQG2+8MWB/YQOViHDhwgVKS0spKSmhtLTU8bzqddXjsWPHyMnJwd3fV/v27QkKCnJ58qzrb7I1McZgjCEoKMjtY9Xzs2fPukwOnhQSElJvIqmZcOr7aej+vPW37okWapMThIsdXYItWaQAF4vIlY2KxMv8KUGUlZU5qrcdO3aMoUOH8tBDD2n1Nh+yWq3VTvA1T/I1H129X99JLDQ0lA4dOnDu3Lk6123fvj3JycluT5g1n9e1rKEn4cZ8TnP3WZUcGspVa6tK1YkvPj6eioqKaj+VlZVUVFRQXl7ueO7up7Kyst713O3PYmleOZygoKBGJaaGrHf8+HHWrl1LZeW/L+w0JUl4JEEYYzpRfdRTsYjU/t/0IX9IEOfPn3dUbyssLGTkyJE89NBDAVW9zVejTKxWK2VlZW5P7DVP8q5O/PX9vrdv357w8HA6dOhAeHh4tefOjzWXVf1U9SM15ISn19Wr8+c+CKvV2uzE5LxdYxKYcwJoiMYes2YlCGPM/cCzQClQtbKISP+Gh9wyfJkgzp07x9tvvx3w1dua06S1Wq0Nujzj7gRfVlZW7wk+LCyszhO8qxO783JP9gf58wnPX2mnfm1Wq9WRNKoeX3jhBc6cOeN2m27duvGnP/2pQftvboLYD1wtIoUN+jQf8kWCqKre9s4773Du3LmArt6WnZ3N3//+dyoq/n1vZEhICNdffz2dO3d2nOyrvuXXPOGXldVfPqQh397dPYaFhfldh7+e8BpP74OonydbqM1NEBnAVBEpadCn+VBLJoia1dsmTJjAgw8+2Cqrt4kIpaWlnDlzhrNnz3LmzJlaPydOnODkyZN17scY06ATu7sEEBYWFpAd93rCU97gqRZqcxPEcOB/gK3AharlIvLLerZbAtwMHBeRwfZlHwBVkXcBTovIMBfb5gLnAAtQ6S74mloiQRw7dsxRva28vNyvq7dZrVbOnTvn9qTvvNy5VVAlNDSUzp0707lzZ44cOeLy20qVrl278l//9V8BeYJXyl95exRTQ6baeBX4GPgOaMw4tDeAhcBbVQtE5HanoF4C3F9Eg+v96bJWXl4er732GsuXL/d59baKigqXJ/2ay9yNpOnQoQOdO3emU6dO9O/f35EEqn46depE586dCQ8Pd3Ss19ekvffeezU5KNXCEhMTmT17ttdaqA1pQXwhItc0aefGxAJrqloQTssNtnsqbhCR/S62ywVGNDZBeKMFUVW9bdWqVRhjvFa9TUQoKyur84Rf9VNSUvtqnzGGjh07ujzR11zWrl27JsWona5KBZ7mtiA2GWNmAaupfomp7gvSdbsWOOYqOVTtHlhvjBHgVRFZ5G5H9thmAfTt27dJwbi6RpyTk8PLL7/sqN42ffp07rvvvkZXb7NarRQXF9d70j979qzLb+chISGOE3vPnj1JSEioddLv3LkzkZGRXu+grfq2op2uSrUNDWlBHHKxuEHDXOtoQbwM5IjIS2626y0i+caYi4ANwCMi8ll9n9fYFkRycjKVlZUMGTKE4OBgLBYL3333HSUlJVRUVBAeHs706dNdVm+rrKx0e7J3TgTu7hANDw93+w3f+XWHDh387v4J7XRVKnB47E7qJnxwLDUShDEmBDgKXCEieQ3YxzPYbsp7sb51G5sgRo0a5UgOVSwWC/v27SMlJYVrrrmGyspKl9/8z58/7ypWIiMj6z3pd+7cucmXeZRSypOadYnJPh/TQ9jqPwi2qb9fEZH6B7W7Nh7Y6y45GGMigCAROWd/PgHbjXoelZ2dXSs5AAQHBzNw4EAOHz7M4cOHAdtlnqqT/EUXXUR8fLzLyzwdO3b0u3H4SinVVA3pg3gL25DTv9tfTwfeBn5a10bGmKXYJvqLMsbkAfNF5HXgDmBpjXV7A4tFZBLQE1hhv6wSArwnIhkN/Qc11BtvvFHnybxjx4785je/oVOnTkRERPjdZR6llPK2hvRBfCMiQ+tb5g8ac4kpOzubP//5zy6ThMVi4fHHH9fr60qpgFfXJaaGDFz/2hhzldPORgGbPRWcryQmJvLdd9/VmqWxqqNak4NSqq1rSIIYBXxhjMm135+wBbjOGPOdMeZbr0bnZSEhIdWSRFVyCAlpUKlupZQKaA05E6Z4PQof2bzZ1hDSYZtKKVWbV4e5tjR/qAehlFKtSXPvpFZKKeWHMjMzWbZsGYWFhURFRTFt2jRGjx7tsf1rglBKqVYoMzOTxYsXO6boKSwsZPHixQAeSxI6/aZSSrVCVSUHnJWXl7Ns2TKPfYa2IJRSqhUQEQoLC8nOziY7O5vCQteTXRcVFXnsMzVBKKWUH7Jarfzwww+OhLBv3z5HVcfw8HBCQ0NdFvrq3r27x2LQBKGUUn6grKyMAwcOOBJCTk4OpaWlAHTr1o3ExEQSExNJSEigb9++fPHFF9X6IMBWn2XatGkei0kThFJK+cDp06fZt2+fIyHk5uZitVoxxhATE0NycrIjIURFRdWaD66qI3rZsmUUFRXRvXt3j49i0vsglFLKy0SE/Pz8agnh2LFjgK32e1xcnKOFMGDAACIiIlosNr0PQimlWlBlZSUHDx509B3s27ePc+fOARAZGUliYiLjxo0jMTGRSy65xG+n9/HPqJRSqhU5f/68o3Wwb98+Dhw44OhA7tmzJ8OHD3e0EHr16tVqygdoglBKqUaoOdx037595OXlISIEBwcTGxvL+PHjHQmhc+fOvg65yTRBKKVUHeobbjpgwABGjRpFYmIicXFxhIWF+Thiz9EEoZRSThoy3DQhIYHExET69u1LUFDgTkihCUIp1aadOXPGkQyys7M5fPgwFoul2nDTqoTgarhpINMEoZRqMxoy3HTy5MmO4aaRkZE+jti3vJYgjDFLgJuB4yIy2L7sGeAXwAn7ak+KyFoX26YAfwOCgcUissBbcSqlAlfVcFPnEUY1h5vecMMNjuGmoaGhPo7Yv3izBfEGsBB4q8byv4rIi+42MsYEA/8AbgTygG3GmFUissdbgSqlAkPVcNOqhOBuuGlCQgK9e/duU5eLmsJrCUJEPjPGxDZh0yuBHBE5CGCMeR9IBTRBKKUcnIebViWEquGmQUFB1YabJiQk0KVLF1+H3Or4og9itjHmbmA7MEdETtV4Pxo44vQ6DxjVUsEppXynrgpp9Q03jY+PD9jhpr7S0gniZeCPgNgfXwJm1ljHVZvP7YRRxphZwCyAvn37eiZKpVSLc1UhbdGiRWzfvp3S0tJaw02rRha1heGmvtKiCUJEjlU9N8a8BqxxsVoe0MfpdQyQX8c+FwGLwDZZn2ciVUq1JBFh6dKltSqkVVZW8tVXX9GnTx+uueYaR0Joa8NNfaVFE4QxppeIFNhf3grsdrHaNmCAMeYS4ChwB/CfLRSiUqoFlJeXk5uby/79+x0/p07VvNr8b88//3wLRqeqeHOY61JgLBBljMkD5gNjjTHDsF0yygXut6/bG9tw1kkiUmmMmQ2swzbMdYmIfO+tOJVS3ldUVFQtGeTm5lJZWQlAjx49GDhwIN988w3nz5+vtW1UVFRLh6vsvDmKabqLxa+7WTcfmOT0ei1Q6/4IpZT/q9k6yMnJcXQmh4aG0r9/f1JSUhgwYADx8fF07doVqN0HAZ6vkKYaR++kVko1S32tg0svvZT4+HgGDBhAv3793NY+aIkKaapxtKKcUqrBGtI6GDBgQK3WgfJfWlFOKdUknmodqNZJ/zeVUgBUVFRw6NChRvcdqMClCUKpNqqu1kFUVJRjRlNtHbRd+j+uVBugrQPVFJoglApA2jpQnqC/FUq1cto6UN6iCUKpVkZbB6ql6G+OUn7MuXWQk5PD/v37tXWgWowmCKX8iLYOlD/R3y6lvKiuAjgVFRW1ZjR11TqouhFNWweqpWmCUMpL3BXAyczMpKSkRFsHyu/pb6BSXlBaWsq7777rsgDOt99+S2JiorYOlN/TBKFUM1VWVvLDDz9w4MABx09+fj51TYQ5f/78FoxQqabRBKFUI1itVgoKChyJ4ODBgxw+fNhxqahTp07ExcVx9dVXs379es6ePVtrH1oAR7UWmiCUckNEOHnyZLWWwaFDhygtLQUgLCzM0ZEcFxdH//79q9VKvuiii7QAjmrVNEEoZVdcXOxoFVQ9nj59GoDg4GD69evH6NGj6d+/P3FxcfTu3ZugoCC3+9MCOKq10wSh2qQLFy6Qm5tbLSEcO3YMAGMMvXv3ZsiQIY6WQb9+/QgNDW3054wePVoTgmq1NEGogGexWDhy5Ei1lsGRI0ewWq0AdO/enbi4OK6//nri4uK45JJL6NChg4+jVsr3vJYgjDFLgJuB4yIy2L7sz8AtQDlwALhXRE672DYXOAdYgEp35fCUqklEOHbsWLWWQW5urqMfICIigri4OIYPH05cXBxxcXF06dLFt0Er5ae82YJ4A1gIvOW0bAMwT0QqjTHPA/OA37rZ/noRKfRifCoAnDp1ypEIqpLC+fPnAVuHcGxsLOPGjXMkg4suusjRiayUqpvXEoSIfGaMia2xbL3Tyy+B//DW56vAU1JSwqFDh6qNKqqamiIoKIg+ffpw5ZVXEh8fT//+/YmJiSE4ONjHUSvVevmyD2Im8IGb9wRYb4wR4FURWeRuJ8aYWcAsgL59+3o8SOUb5eXl1W4+O3jwIPn5+Y73e/bsyaWXXupoGfTr14/27dv7MGKlAo9PEoQx5imgEnjXzSrJIpJvjLkI2GCM2Ssin7la0Z48FgGMGDHC/a2rym9ZrVby8/PJyclxXC764YcfsFgsAHTp0oW4uDjHENP+/fsTGRnp46iVCnwtniCMMfdg67weJ27mIhCRfPvjcWPMCuBKwGWCUK2LiFBYWFitZXDo0CHKysoACA8Pp3///kyePNlxv0G3bt2030ApH2jRBGGMScHWKX2diJS4WScCCBKRc/bnE4BnWzBM5UFnz56tdfNZ1fQTISEhxMbGMmbMGMf9Br169arz5jOlVMvx5jDXpcBYIMoYkwfMxzZqqT22y0YAX4rIA8aY3sBiEZkE9ARW2N8PAd4TkQxvxakap676BmVlZRw6dKjaqKITJ04AtpvPoqOjGT58uKNl0LdvX53SWik/ZuqacbK1GTFihGzfvt3XYQSsmvUNwDYFRUJCAsXFxeTl5TlmMO3Ro4cjEcTFxREbG0t4eLivQldKuWGM2eHuXjP9+qYapKysjHfeeadWfQOLxcLevXsZOnQoI0eOdFwq6ty5s48iVUp5iiYI5VJZWRn79+9nz549ZGVlceDAAceooppEhCeeeKKFI1RKeZsmCAW4TwjBwcGOUUWffPKJ1jdQqg3RBNFGlZWVsW/fPkdCOHjwYK2EkJSUREJCAmFhYQDExMRofQOl2hBNEG1EUxJCTVrfQKm2RUcxBaj6EsLAgQPrTQhKqcCno5jagLoSQlxcHDfffDNJSUkMGDBAE4JSqkE0QbRSmhCUUt6mCaKVKCsrIzs7m6ysLPbs2cOhQ4c0ISilvEoThJ/ShKCU8jVNEH6iZkI4ePAgVqtVE0IAuHDhAnv27CEpKUlrVjSAHi//oQnCR+pLCFOmTGHgwIF+mRD0D7hxcnNzOXPmDLm5uSQmJvo6HL+nx8t/aIJoIaWlpbU6lVtLQqjJ3/6ARQQRwWq1euS5p/ZV867zgoICCgoKAOjQoYPLGhc1l7XGdZq6XdWxcX5dUFCAMYbY2FiMMQQFBdX52JB1nNcNhDoj3vzCpgnCSwIpIYCt6tvnn3+O830zzn/AiYmJPjsxtwTnE4qrk1HN58HBwXTp0oWSkpJqd563b9+eiIiIajUvXN2L1JBlda1T87Gp+2nJdUJCQrBYLC63OXToUK3tPKG5ScYTiaox27rizS9smiA8pLUkBIvFQmVlJRUVFbUeXS2reqzrRCwi7N27t87PbeiJ1fkPIiQkpFHbNOV5Y9ZriuzsbAoKCggKCsJqtdKtWze/aHX5q5rHq1evXiQkJLj9ktDUx6ZsU5W86tuHNzn/fVRUVFR7r+oLW1BQEGPGjPHI52mCaKK6EkJ8fDxTpkxxdCp74zq9N070xhhCQ0MJCQkhNDSUsLAwQkNDqy07ceIEp06dwhiDiNCjRw9iY2PrPcm2VRUVFfTu3ZvevXuTn59fa7p0VZ2r4+X8u9QaeDtRVT1WVFRw9uxZLly4ANiSR1RUFHFxcR77t2iCaCBvJQRfnehrPoaGhjboZH7y5Mlaf8AREREN/ve2NYMHD3Y8T0hI8GEkrUMgHK+qS4wtoWaLKzg42KNfSNt8gsjMzGTVqlUMGzaMr7/+mtTUVEaPHk1paWmt+xDqSgjOJ/rS0lKvnuidT+quTvjBwcFe+9YeCH/ASgUKb7dQ2/RkfVUlNEeOHElSUhJ79uxh27ZtXHzxxZw/f5527doRHh5Onz59iI6OpkePHnTs2BGr1dqsE72rb+8tfaJXSinQyfrcunDhAjNmzHC8HjRoEIMGDXK7fllZGRaLxa++0SullLd4LUEYY5YANwPHRWSwfVk34AMgFsgFponIKRfbpgB/A4KBxSKywBsxvvfee1x11VXExcU5ruGdPHmSnJwcZs6cqSd6pVSb5s1hAW8AKTWWzQU2isgAYKP9dTXGmGDgH8BNQBIw3RiT5I0AIyIiqKiowBhDZWUlxhiOHTtGfn4+F110EV27dqVjx46EhYUREhKiyUEp1aZ4LUGIyGfAyRqLU4E37c/fBH7iYtMrgRwROSgi5cD79u08btq0aURERLBnzx5WrlzJnj17iIyM1BKaSimFd1sQrvQUkQIA++NFLtaJBo44vc6zL3PJGDPLGLPdGLP9xIkTjQpm9OjRxMXFkZ2dzalTp8jOziYuLk5LaCqlFP7ZSe3qOo7boVYisghYBLZRTI39sNGjR2tCUEopF1q6BXHMGNMLwP543MU6eUAfp9cxQH4LxKaUUspJSyeIVcA99uf3AOku1tkGDDDGXGKMaQfcYd9OKaVUC/JagjDGLAW2AInGmDxjzM+BBcCNxpj9wI321xhjehtj1gKISCUwG1gHZAHLROR7b8WplFLKNa/1QYjIdDdvjXOxbj4wyen1WmCtl0JTSinVAK1jekSllFItLqDmYjLGnAAON3HzKKDQg+EEOj1ejaPHq3H0eDVOc45XPxHp4eqNgEoQzWGM2e5uwipVmx6vxtHj1Th6vBrHW8dLLzEppZRySROEUkoplzRB/NsiXwfQyujxahw9Xo2jx6txvHK8tA9CKaWUS9qCUEop5ZImCKWUUi4FbIIwxvQxxmwyxmQZY743xjxqX97NGLPBGLPf/tjVvry7ff1iY8xCp/10NMbscvopNMb8Xx/9s7ymCcfrSqdj8o0x5lanfV1hjPnOGJNjjPlvE4CVlhp7vJy262v/HXvMaZker9q/X7HGmFKn37FXnPalx8vF75cx5jJjzBb7+t8ZY8Lsy5t+vEQkIH+AXsDl9ucdgX3YKtS9AMy1L58LPG9/HgGMBh4AFtax3x3AGF//+/zgeHUAQpy2Pe70+ivgamxTt38I3OTrf5+vj5fTdsuB/wc85rRMj1ft369YYLebfenxqn28QoBvgaH2192B4OYer4BtQYhIgYjstD8/h23iv2jcVLUTkfMikgmUudunMWYAtiJHn3svct9owvEqEdvEigBh2Gt22Kdx7yQiW8T22/kWrisHtmqNPV4AxpifAAeB752W6fGqv8qkgx4vt8drAvCtiHxj36ZIRCzNPV4BmyCcGWNigeHAVhpW1c6d6cAH9gMdsBp6vIwxo4wx3wPfAQ/YE0Y0tpoeVeqsCBgIGnK8jDERwG+BP9TYXI+X+7/HS4wxXxtjPjXGXGtfpsfL9fFKAMQYs84Ys9MY84R9ebOOlz9WlPMoY0wktmb9r0TkbDMvV94B3OWRwPxUY46XiGwFBhljBgJvGmM+pJEVAVu7RhyvPwB/FZHiGuvo8XKtAOgrIkXGmCuAlcaYQejxcrdqCLZL5COBEmCjMWYHcNbFug0+XgGdIIwxodgO7rsikmZffMwY00tECoz7qnau9jUU2zX2HV4K1+eaerxEJMsYcx4YjO0bSozT2wFbEbCRx2sU8B/GmBeALoDVGFNm316PV43jJSIXgAv25zuMMQewfUvW3y/Xv195wKciUmjfdi1wOfAOzTheAXuJyd5T/zqQJSJ/cXqrIVXtXJkOLPVchP6lscfL2Cr+hdif9wMSgVx7s/ecMeYq+z7vpuHHuNVo7PESkWtFJFZEYoH/C/yXiCzU4+X296uHMSbY/rw/MAA4qMfL7flrHXCZMaaD/e/yOmBPs49XS/fOt9QPtuaWYOvZ32X/mYStd38jsN/+2M1pm1zgJFCMLSMnOb13ELjU1/8ufzle2C61fW9fbyfwE6d9jQB2AweAhdjv2A+kn6b8fjlt+wzVRzHp8ar9+3Wb/ffrG/vv1y16vOo9f/3Mfsx2Ay944njpVBtKKaVcCthLTEoppZpHE4RSSimXNEEopZRySROEUkoplzRBKKWUckkThFJKKZc0QSillHLp/wOJFNnHdjJ5dQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmoysf       = np.zeros((2,4))\n",
    "pmoys0f      = np.zeros((2,4))\n",
    "pmoysf[0,:]  = yy_eu\n",
    "pmoysf[1,:]  = yy_us\n",
    "pmoys0f[0,:] = yy_eu0\n",
    "pmoys0f[1,:] = yy_us0\n",
    "\n",
    "\n",
    "xx      = np.zeros((2,5))\n",
    "xx[0,0] = PMoY10[0]\n",
    "xx[1,0] = PMoY10[-1]\n",
    "zz      = np.zeros((2,5))\n",
    "zz[0,0] = PMoY10[0]\n",
    "zz[1,0] = PMoY10[-1]\n",
    "for jj in range(2):\n",
    "    for ii in range(horizon):\n",
    "        xx[jj,ii+1] = pmoysf[jj,ii] \n",
    "        zz[jj,ii+1] = pmoys0f[jj,ii]\n",
    "    \n",
    "lyhorizon = ('2017','2030','2040','2050','2060')\n",
    "\n",
    "plt.plot(np.arange(len(lyhorizon)),list(xx[0,:]),'o-',label=r'EU, $\\Delta p \\neq 0$',color = '0.35')\n",
    "plt.plot(np.arange(len(lyhorizon)),list(xx[1,:]),'s-',label=r'US, $\\Delta p \\neq 0$',color = '0.15')\n",
    "plt.plot(np.arange(len(lyhorizon)),list(zz[0,:]),'*-',label=r'EU, $\\Delta p = 0$',color = '0.75')\n",
    "plt.plot(np.arange(len(lyhorizon)),list(zz[1,:]),'D-',label=r'US, $\\Delta p = 0$',color = '0.4')\n",
    "plt.xticks(np.arange(len(lyhorizon)), labels=lyhorizon)\n",
    "plt.legend()\n",
    "plt.ylabel('pm/y (%)')\n",
    "plt.savefig('../figures/fig_a2_forecast.eps')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f543968e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
